#
# kppm.R
#
# kluster/kox point process models
#
# $Revision: 1.158 $ $Date: 2020/11/16 01:32:06 $
#

kppm <- function(X, ...) {
  UseMethod("kppm")
}


kppm.formula <-
  function(X, clusters = c("Thomas","MatClust","Cauchy","VarGamma","LGCP"),
           ..., data=NULL) {
  ## remember call
  callstring <- short.deparse(sys.call())
  cl <- match.call()

  ########### INTERPRET FORMULA ##############################
  
  if(!inherits(X, "formula"))
    stop(paste("Argument 'X' should be a formula"))
  formula <- X
  
  if(spatstat.options("expand.polynom"))
    formula <- expand.polynom(formula)

  ## check formula has LHS and RHS. Extract them
  if(length(formula) < 3)
    stop(paste("Formula must have a left hand side"))
  Yexpr <- formula[[2L]]
  trend <- formula[c(1L,3L)]
  
  ## FIT #######################################
  thecall <- call("kppm", X=Yexpr, trend=trend,
                  data=data, clusters=clusters)
  ncall <- length(thecall)
  argh <- list(...)
  nargh <- length(argh)
  if(nargh > 0) {
    thecall[ncall + 1:nargh] <- argh
    names(thecall)[ncall + 1:nargh] <- names(argh)
  }
##  result <- eval(thecall, 
##                 envir=if(!is.null(data)) data else parent.frame(),
##                 enclos=if(!is.null(data)) parent.frame() else baseenv())
  callenv <- list2env(as.list(data), parent=parent.frame())
  result <- eval(thecall, envir=callenv, enclos=baseenv())

  result$call <- cl
  result$callframe <- parent.frame()
  if(!("callstring" %in% names(list(...))))
    result$callstring <- callstring
  
  return(result)
}

kppm.ppp <- kppm.quad <-
  function(X, trend = ~1,
           clusters = c("Thomas","MatClust","Cauchy","VarGamma","LGCP"),
           data=NULL,
           ...,
           covariates = data,
           subset, 
           method = c("mincon", "clik2", "palm"),
           improve.type = c("none", "clik1", "wclik1", "quasi"),
           improve.args = list(),
           weightfun=NULL,
           control=list(),
           algorithm="Nelder-Mead",
           statistic="K",
           statargs=list(),
           rmax = NULL,
           covfunargs=NULL,
           use.gam=FALSE,
           nd=NULL, eps=NULL) {
  cl <- match.call()
  callstring <- paste(short.deparse(sys.call()), collapse="")
  Xname <- short.deparse(substitute(X))
  clusters <- match.arg(clusters)
  improve.type <- match.arg(improve.type)
  method <- match.arg(method)
  if(method == "mincon")
    statistic <- pickoption("summary statistic", statistic,
                            c(K="K", g="pcf", pcf="pcf"))
  ClusterArgs <- list(method = method,
                      improve.type = improve.type,
                      improve.args = improve.args,
                      weightfun=weightfun,
                      control=control,
                      algorithm=algorithm,
                      statistic=statistic,
                      statargs=statargs,
                      rmax = rmax)
  Xenv <- list2env(as.list(covariates), parent=parent.frame())
  X <- eval(substitute(X), envir=Xenv, enclos=baseenv())
  isquad <- is.quad(X)
  if(!is.ppp(X) && !isquad)
    stop("X should be a point pattern (ppp) or quadrature scheme (quad)")
  if(is.marked(X))
    stop("Sorry, cannot handle marked point patterns")
  if(!missing(subset)) {
    W <- eval(subset, covariates, parent.frame())
    if(!is.null(W)) {
      if(is.im(W)) {
        W <- solutionset(W)
      } else if(!is.owin(W)) {
        stop("Argument 'subset' should yield a window or logical image",
             call.=FALSE)
      }
      X <- X[W]
    }
  }
  po <- ppm(Q=X, trend=trend, covariates=covariates,
            forcefit=TRUE, rename.intercept=FALSE,
            covfunargs=covfunargs, use.gam=use.gam, nd=nd, eps=eps)
  XX <- if(isquad) X$data else X
  # set default weight function
  if(is.null(weightfun) && method != "mincon") {
    RmaxW <- (rmax %orifnull% rmax.rule("K", Window(XX), intensity(XX))) / 2
    weightfun <- function(d, rr=RmaxW) { as.integer(d <= rr) }
    formals(weightfun)[[2]] <- RmaxW
    attr(weightfun, "selfprint") <- paste0("Indicator(distance <= ", RmaxW, ")")
  }
  # fit
  out <- switch(method,
         mincon = kppmMinCon(X=XX, Xname=Xname, po=po, clusters=clusters,
                             control=control, statistic=statistic,
                             statargs=statargs, rmax=rmax,
                             algorithm=algorithm, ...),
         clik2   = kppmComLik(X=XX, Xname=Xname, po=po, clusters=clusters,
                             control=control, weightfun=weightfun, 
                             rmax=rmax, algorithm=algorithm, ...),
         palm   = kppmPalmLik(X=XX, Xname=Xname, po=po, clusters=clusters,
                             control=control, weightfun=weightfun, 
                             rmax=rmax, algorithm=algorithm, ...))
  ##
  h <- attr(out, "h")
  out <- append(out, list(ClusterArgs=ClusterArgs,
                          call=cl,
                          callframe=parent.frame(),
                          callstring=callstring))
  # Detect DPPs
  DPP <- list(...)$DPP
  class(out) <- c(ifelse(is.null(DPP), "kppm", "dppm"), class(out))

  # Update intensity estimate with improve.kppm if necessary:
  if(improve.type != "none")
    out <- do.call(improve.kppm,
                   append(list(object = out, type = improve.type),
                          improve.args))

  attr(out, "h") <- h
  return(out)
}

kppmMinCon <- function(X, Xname, po, clusters, control, statistic, statargs,
                       algorithm="Nelder-Mead", DPP=NULL, ...) {
  # Minimum contrast fit
  stationary <- is.stationary(po)
  # compute intensity
  if(stationary) {
    lambda <- summary(po)$trend$value
  } else {
    # compute intensity at high resolution if available
    w <- as.owin(po, from="covariates")
    if(!is.mask(w)) w <- NULL
    lambda <- predict(po, locations=w)
  }
  # Detect DPP model and change clusters and intensity correspondingly
  if(!is.null(DPP)){
    tmp <- dppmFixIntensity(DPP, lambda, po)
    clusters <- tmp$clusters
    lambda <- tmp$lambda
    po <- tmp$po
  }
  mcfit <- clusterfit(X, clusters, lambda = lambda,
                      dataname = Xname, control = control,
                      statistic = statistic, statargs = statargs,
                      algorithm=algorithm, ...)
  fitinfo <- attr(mcfit, "info")
  attr(mcfit, "info") <- NULL
  # all info that depends on the fitting method:
  Fit <- list(method    = "mincon",
              statistic = statistic,
              Stat      = fitinfo$Stat,
              StatFun   = fitinfo$StatFun,
              StatName  = fitinfo$StatName,
              FitFun    = fitinfo$FitFun,
              statargs  = statargs,
              mcfit     = mcfit,
              maxlogcl  = NULL)
  # results
  if(!is.null(DPP)){
    clusters <- update(clusters, as.list(mcfit$par))
    out <- list(Xname      = Xname,
                X          = X,
                stationary = stationary,
                fitted     = clusters,
                po         = po,
                Fit        = Fit)
  } else{
    out <- list(Xname      = Xname,
                X          = X,
                stationary = stationary,
                clusters   = clusters,
                modelname  = fitinfo$modelname,
                isPCP      = fitinfo$isPCP,
                po         = po,
                lambda     = lambda,
                mu         = mcfit$mu,
                par        = mcfit$par,
                par.canon  = mcfit$par.canon,
                clustpar   = mcfit$clustpar,
                clustargs  = mcfit$clustargs,
                modelpar   = mcfit$modelpar,
                covmodel   = mcfit$covmodel,
                Fit        = Fit)
  }
  attr(out, "h") <- attr(mcfit, "h")
  return(out)
}

clusterfit <- function(X, clusters, lambda = NULL, startpar = NULL,
                       ...,
                       q=1/4, p=2, rmin=NULL, rmax=NULL, 
                       ctrl=list(q=q, p=p, rmin=rmin, rmax=rmax),
                       statistic = NULL, statargs = NULL,
                       algorithm="Nelder-Mead", verbose=FALSE,
                       pint=NULL){
  if(verbose) splat("Fitting cluster model")
  ## If possible get dataname from dots
  dataname <- list(...)$dataname
  ## Cluster info:
  info <- spatstatClusterModelInfo(clusters)
  if(verbose) splat("Retrieved cluster model information")
  ## Determine model type
  isPCP <- info$isPCP
  isDPP <- inherits(clusters, "detpointprocfamily")

  ## resolve algorithm parameters
  default.ctrl <- list(q=if(isDPP) 1/2 else 1/4,
                       p=2,
                       rmin=NULL,
                       rmax=NULL)
  given.ctrl <- if(missing(ctrl)) list() else ctrl[names(default.ctrl)]
  given.args <- c(if(missing(q)) NULL else list(q=q),
                  if(missing(p)) NULL else list(p=p),
                  if(missing(rmin)) NULL else list(rmin=rmin),
                  if(missing(rmax)) NULL else list(rmax=rmax))
  ctrl <- resolve.defaults(given.args, given.ctrl, default.ctrl)
  if(verbose) {
    splat("Algorithm parameters:")
    print(ctrl)
  }
  ##
  if(inherits(X, "ppp")){
      if(verbose) 
        splat("Using point pattern data")
      if(is.null(dataname))
         dataname <- getdataname(short.deparse(substitute(X), 20), ...)
      if(is.null(statistic))
          statistic <- "K"
      # Startpar:
      if(is.null(startpar))
          startpar <- info$selfstart(X)
      stationary <- is.null(lambda) || (is.numeric(lambda) && length(lambda)==1)
      if(verbose) {
        splat("Starting parameters:")
        print(startpar)
        cat("Calculating summary function...")
      }
      # compute summary function
      if(stationary) {
          if(is.null(lambda)) lambda <- intensity(X)
          StatFun <- if(statistic == "K") "Kest" else "pcf"
          StatName <-
              if(statistic == "K") "K-function" else "pair correlation function"
          Stat <- do.call(StatFun,
                          resolve.defaults(list(X=quote(X)),
                                           statargs,
                                           list(correction="best")))
      } else {
          StatFun <- if(statistic == "K") "Kinhom" else "pcfinhom"
          StatName <- if(statistic == "K") "inhomogeneous K-function" else
          "inhomogeneous pair correlation function"
          Stat <- do.call(StatFun,
                          resolve.defaults(list(X=quote(X), lambda=lambda),
                                           statargs,
                                           list(correction="best")))
      }
      if(verbose) splat("Done.")
  } else if(inherits(X, "fv")){
      if(verbose) 
        splat("Using the given summary function")
      Stat <- X
      ## Get statistic type
      stattype <- attr(Stat, "fname")
      StatFun <- paste0(stattype)
      StatName <- NULL
      if(is.null(statistic)){
          if(is.null(stattype) || !is.element(stattype[1L], c("K", "pcf")))
              stop("Cannot infer the type of summary statistic from argument ",
                   sQuote("X"), " please specify this via argument ",
                   sQuote("statistic"))
          statistic <- stattype[1L]
      }
      if(stattype[1L]!=statistic)
          stop("Statistic inferred from ", sQuote("X"),
               " not equal to supplied argument ",
               sQuote("statistic"))
      # Startpar:
      if(is.null(startpar)){
          if(isDPP)
              stop("No rule for starting parameters in this case. Please set ",
                   sQuote("startpar"), " explicitly.")
          startpar <- info$checkpar(startpar, old=FALSE)
          startpar[["scale"]] <- mean(range(Stat[[fvnames(Stat, ".x")]]))
      }
  } else{
      stop("Unrecognised format for argument X")
  }
  
  ## avoid using g(0) as it may be infinite
  if(statistic=="pcf"){
      if(verbose) splat("Checking g(0)")
      argu <- fvnames(Stat, ".x")
      rvals <- Stat[[argu]]
      if(rvals[1L] == 0 && (is.null(rmin) || rmin == 0)) {
          if(verbose) splat("Ignoring g(0)")
          rmin <- rvals[2L]
      }
  }

  ## DPP resolving algorithm and checking startpar
  changealgorithm <- length(startpar)==1 && algorithm=="Nelder-Mead"
  if(isDPP){
    if(verbose) splat("Invoking dppmFixAlgorithm")
    alg <- dppmFixAlgorithm(algorithm, changealgorithm, clusters, startpar)
    algorithm <- alg$algorithm
  }

  dots <- info$resolvedots(...)
  #' determine initial values of parameters
  startpar <- info$checkpar(startpar)
  #' code to compute the theoretical summary function of the model
  theoret <- info[[statistic]]
  #' explanatory text
  desc <- paste("minimum contrast fit of", info$descname)

  #' ............ experimental .........................
  strict <- !isFALSE(spatstat.options("kppm.strict"))
  sargs <- if(strict) list() else list(strict=FALSE)
  
  #' ............ experimental .........................
  do.adjust <- spatstat.options("kppm.adjusted")
  if(do.adjust) {
    if(verbose) splat("Applying kppm adjustment")
    W <- Window(X)
    adjdata <- list(paircorr = info[["pcf"]],
                    pairWcdf = distcdf(W),
		    tohuman  = NULL)
    adjfun <- function(theo, par, auxdata, ...) {
      with(auxdata, {
        if(!is.null(tohuman))
	  par <- tohuman(par)
        a <- as.numeric(stieltjes(paircorr, pairWcdf, par=par, ...))
	return(theo/a)
      })
    }
    adjustment <- list(fun=adjfun, auxdata=adjdata)
  } else adjustment <- NULL
    
  #' ............ experimental .........................
  usecanonical <- spatstat.options("kppm.canonical")
  if(usecanonical) {
    if(verbose) splat("Converting to canonical parameters")
     tocanonical <- info$tocanonical
     tohuman <- info$tohuman
     if(is.null(tocanonical) || is.null(tohuman)) {
       warning("Canonical parameters are not yet supported for this model")
       usecanonical <- FALSE
     } 

  }
  startpar.human <- startpar
  if(usecanonical) {
    htheo <- theoret
    startpar <- tocanonical(startpar)
    theoret <- function(par, ...) { htheo(tohuman(par), ...) }
    if(do.adjust)
      adjustment$auxdata$tohuman <- tohuman
  }
  #' ............ experimental .........................
  whiu <- pint$whiu
  if(is.function(whiu) && usecanonical) {
    whiu.human <- whiu
    whiu <- function(par, ...) { whiu.human(tohuman(par), ...) }
    pint$whiu <- whiu
  }
  #' ...................................................
  
  #'
  mcargs <- resolve.defaults(list(observed=Stat,
                                  theoretical=theoret,
                                  startpar=startpar,
                                  ctrl=ctrl,
                                  method=algorithm,
                                  fvlab=list(label="%s[fit](r)",
                                      desc=desc),
                                  explain=list(dataname=dataname,
                                      fname=statistic,
                                      modelname=info$modelname),
                                  margs=dots$margs,
                                  model=dots$model,
                                  funaux=info$funaux,
				  adjustment=adjustment,
                                  pint=pint),
                             list(...),
                             sargs)

  if(isDPP && algorithm=="Brent" && changealgorithm)
    mcargs <- resolve.defaults(mcargs, list(lower=alg$lower, upper=alg$upper))

  ## .............. FIT .......................
  if(verbose) splat("Starting minimum contrast fit")
  mcfit <- do.call(mincontrast, mcargs)
  if(verbose) splat("Returned from minimum contrast fit")
  ## ..........................................
  
  ## extract fitted parameters and reshape
  if(!usecanonical) {
    optpar.canon <- NULL
    optpar.human <- mcfit$par
    names(optpar.human) <- names(startpar.human)
  } else {
    optpar.canon <- mcfit$par
    names(optpar.canon) <- names(startpar)
    optpar.human <- tohuman(optpar.canon)
    names(optpar.human) <- names(startpar.human)
  }
  mcfit$par       <- optpar.human
  mcfit$par.canon <- optpar.canon
  # Return results for DPPs
  if(isDPP){
    extra <- list(Stat      = Stat,
                  StatFun   = StatFun,
                  StatName  = StatName,
                  modelname  = info$modelabbrev,
                  lambda     = lambda)
    attr(mcfit, "info") <- extra
    if(verbose) splat("Returning from clusterfit (DPP case)")
    return(mcfit)
  }
  ## Extra stuff for ordinary cluster/lgcp models
  ## imbue with meaning
  ## infer model parameters
  mcfit$modelpar <- info$interpret(optpar.human, lambda)
  mcfit$internal <- list(model=ifelse(isPCP, clusters, "lgcp"))
  mcfit$covmodel <- dots$covmodel
  
  if(isPCP) {
    # Poisson cluster process: extract parent intensity kappa
    kappa <- mcfit$par[["kappa"]]
    # mu = mean cluster size
    mu <- lambda/kappa
  } else {
    # LGCP: extract variance parameter sigma2
    sigma2 <- mcfit$par[["sigma2"]]
    # mu = mean of log intensity 
    mu <- log(lambda) - sigma2/2
  }
  ## Parameter values (new format)
  mcfit$mu <- mu
  mcfit$clustpar <- info$checkpar(mcfit$par, old=FALSE, strict=strict)
  mcfit$clustargs <- info$checkclustargs(dots$margs, old=FALSE, strict=strict)

  ## The old fit fun that would have been used (DO WE NEED THIS?)
  FitFun <- paste0(tolower(clusters), ".est", statistic)

  extra <- list(FitFun    = FitFun,
                Stat      = Stat,
                StatFun   = StatFun,
                StatName  = StatName,
                modelname  = info$modelabbrev,
                isPCP      = isPCP,
                lambda     = lambda)
  attr(mcfit, "info") <- extra
  if(verbose) splat("Returning from clusterfit")
  return(mcfit)
}


kppmComLik <- function(X, Xname, po, clusters, control, weightfun, rmax,
                       algorithm="Nelder-Mead", DPP=NULL, ..., pint=NULL) {
  W <- as.owin(X)
  if(is.null(rmax))
    rmax <- rmax.rule("K", W, intensity(X))
  # identify pairs of points that contribute
  cl <- closepairs(X, rmax, what="ijd")
#  I <- cl$i
#  J <- cl$j
  dIJ <- cl$d
  # compute weights for pairs of points
  if(is.function(weightfun)) {
    wIJ <- weightfun(dIJ)
    sumweight <- safevalue(sum(wIJ))
  } else {
    npairs <- length(dIJ)
    wIJ <- rep.int(1, npairs)
    sumweight <- npairs
  }
  # convert window to mask, saving other arguments for later
  dcm <- do.call.matched(as.mask,
                         append(list(w=W), list(...)),
                         sieve=TRUE)
  M         <- dcm$result
  otherargs <- dcm$otherargs

  ## Detect DPP usage
  isDPP <- inherits(clusters, "detpointprocfamily")

  # compute intensity at pairs of data points
  # and c.d.f. of interpoint distance in window
  if(stationary <- is.stationary(po)) {
    # stationary unmarked Poisson process
    lambda <- intensity(X)
#    lambdaIJ <- lambda^2
    # compute cdf of distance between two uniform random points in W
    g <- distcdf(W)
    # scaling constant is (area * intensity)^2
    gscale <- npoints(X)^2  
  } else {
    # compute fitted intensity at data points and in window
#    lambdaX <- fitted(po, dataonly=TRUE)
    lambda <- lambdaM <- predict(po, locations=M)
    # lambda(x_i) * lambda(x_j)
#    lambdaIJ <- lambdaX[I] * lambdaX[J]
    # compute cdf of distance between two random points in W
    # with density proportional to intensity function
    g <- distcdf(M, dW=lambdaM)
    # scaling constant is (integral of intensity)^2
    gscale <- safevalue(integral.im(lambdaM)^2, default=npoints(X)^2)
  }

  # Detect DPP model and change clusters and intensity correspondingly
  isDPP <- !is.null(DPP)
  if(isDPP){
    tmp <- dppmFixIntensity(DPP, lambda, po)
    clusters <- tmp$clusters
    lambda <- tmp$lambda
    po <- tmp$po
  }

  # trim 'g' to [0, rmax] 
  g <- g[with(g, .x) <= rmax,]
  # get pair correlation function (etc) for model
  info <- spatstatClusterModelInfo(clusters)
  pcfun      <- info$pcf
  funaux     <- info$funaux
  selfstart  <- info$selfstart
  isPCP      <- info$isPCP
  parhandler <- info$parhandler
  modelname  <- info$modelname
  # Assemble information required for computing pair correlation
  pcfunargs <- list(funaux=funaux)
  if(is.function(parhandler)) {
    # Additional parameters of cluster model are required.
    # These may be given as individual arguments,
    # or in a list called 'covmodel'
    clustargs <- if("covmodel" %in% names(otherargs))
                 otherargs[["covmodel"]] else otherargs
    clargs <- do.call(parhandler, clustargs)
    pcfunargs <- append(clargs, pcfunargs)
  } else clargs <- NULL
  # determine starting parameter values
  startpar <- selfstart(X)
  #' ............ experimental .........................
  strict <- !isFALSE(spatstat.options("kppm.strict"))
  if(!strict) pcfunargs <- append(pcfunargs, list(strict=FALSE))
  #' ............ experimental .........................
  usecanonical <- spatstat.options("kppm.canonical")
  if(usecanonical) {
     tocanonical <- info$tocanonical
     tohuman <- info$tohuman
     if(is.null(tocanonical) || is.null(tohuman)) {
       warning("Canonical parameters are not yet supported for this model")
       usecanonical <- FALSE
     }
  }
  startpar.human <- startpar
  if(usecanonical) {
    pcftheo <- pcfun
    startpar <- tocanonical(startpar)
    pcfun <- function(par, ...) { pcftheo(tohuman(par), ...) }
  }
  #' ............ experimental/debugger .........................
  whiu <- pint$whiu
  if(is.function(whiu) && usecanonical) {
    whiu.human <- whiu
    whiu <- function(par, ...) { whiu.human(tohuman(par), ...) }
  }
  TRACE <- isTRUE(pint$trace)
  if(SAVE <- isTRUE(pint$save)) {
    saveplace <- new.env()
    assign("h", NULL, envir=saveplace)
  } else saveplace <- NULL
  # .....................................................
  # create local function to evaluate pair correlation
  #  (with additional parameters 'pcfunargs' in its environment)
  paco <- function(d, par) {
    do.call(pcfun, append(list(par=par, rvals=d), pcfunargs))
  }
  #' ..........  define objective function ......................
  if(!is.function(weightfun)) {
    # pack up necessary information
    objargs <- list(dIJ=dIJ, sumweight=sumweight, g=g, gscale=gscale, 
                    envir=environment(paco),
                    whiu=NULL,   # updated below
                    TRACE=FALSE, # updated below
                    saveplace=NULL, # updated below
                    BIGVALUE=1, # updated below
                    SMALLVALUE=.Machine$double.eps)
    # define objective function (with 'paco' in its environment)
    # This is the log composite likelihood minus the constant term 
    #       sum(log(lambdaIJ)) - npairs * log(gscale)
    obj <- function(par, objargs) {
      with(objargs, {
        logprod <- sum(log(safevalue(paco(dIJ, par))))
        integ <- unlist(stieltjes(paco, g, par=par))
        integ <- pmax(SMALLVALUE, integ)
        logcl <- 2*(logprod - sumweight * log(integ))
        if(is.function(whiu)) logcl <- whiu(par, logcl, 1)
        logcl <- safevalue(logcl, default=-BIGVALUE)
        ## debugger
        if(isTRUE(TRACE)) {
          cat("Parameters:", fill=TRUE)
          print(par)
          splat("log composite likelihood:", logcl)
        }
        if(is.environment(saveplace)) {
          h <- get("h", envir=saveplace)
          hplus <- as.data.frame(append(par, list(logcl=logcl)))
          h <- rbind(h, hplus)
          assign("h", h, envir=saveplace)
        }
        return(logcl)
      },
      enclos=objargs$envir)
    }
    ## determine a suitable large number to replace Inf
    objargs$BIGVALUE <- bigvaluerule(obj, objargs, startpar)
    objargs$whiu <- whiu
    objargs$saveplace <- saveplace
    objargs$TRACE <- TRACE
  } else {
    # create local function to evaluate  pair correlation(d) * weight(d)
    #  (with additional parameters 'pcfunargs', 'weightfun' in its environment)
    force(weightfun)
    wpaco <- function(d, par) {
      y <- do.call(pcfun, append(list(par=par, rvals=d), pcfunargs))
      w <- weightfun(d)
      return(y * w)
    }
    # pack up necessary information
    objargs <- list(dIJ=dIJ, wIJ=wIJ, sumweight=sumweight, g=g, gscale=gscale, 
                    envir=environment(wpaco),
                    whiu=NULL,   # updated below
                    TRACE=FALSE, # updated below
                    saveplace=NULL, # updated below
                    BIGVALUE=1, # updated below
                    SMALLVALUE=.Machine$double.eps)
    # define objective function (with 'paco', 'wpaco' in its environment)
    # This is the log composite likelihood minus the constant term 
    #       sum(wIJ * log(lambdaIJ)) - sumweight * log(gscale)
    obj <- function(par, objargs) {
      with(objargs,
      {
        integ <- unlist(stieltjes(wpaco, g, par=par))
        integ <- pmax(SMALLVALUE, integ)
        logcl <- safevalue(
          2*(sum(wIJ * log(safevalue(paco(dIJ, par))))
            - sumweight * log(integ)),
          default=-BIGVALUE)
        if(is.function(whiu)) logcl <- whiu(par, logcl, 1)
        ## debugger
        if(isTRUE(TRACE)) {
          cat("Parameters:", fill=TRUE)
          print(par)
          splat("log composite likelihood:", logcl)
        }
        if(is.environment(saveplace)) {
          h <- get("h", envir=saveplace)
          hplus <- as.data.frame(append(par, list(logcl=logcl)))
          h <- rbind(h, hplus)
          assign("h", h, envir=saveplace)
        }
        return(logcl)
      },
      enclos=objargs$envir)
    }
    ## determine a suitable large number to replace Inf
    objargs$BIGVALUE <- bigvaluerule(obj, objargs, startpar)
    objargs$whiu <- whiu
    objargs$saveplace <- saveplace
    objargs$TRACE <- TRACE
  }
  #' .........................................................
  ## arguments for optimization
  ctrl <- resolve.defaults(list(fnscale=-1), control, list(trace=0))
  optargs <- list(par=startpar, fn=obj, objargs=objargs, control=ctrl, method=algorithm)
  ## DPP resolving algorithm and checking startpar
  changealgorithm <- length(startpar)==1 && algorithm=="Nelder-Mead"
  if(isDPP){
    alg <- dppmFixAlgorithm(algorithm, changealgorithm, clusters,
                            startpar.human)
    algorithm <- optargs$method <- alg$algorithm
    if(algorithm=="Brent" && changealgorithm){
      optargs$lower <- alg$lower
      optargs$upper <- alg$upper
    }
  }
  ## ..........   optimize it ..............................
  opt <- do.call(optim, optargs)

  ## raise warning/error if something went wrong
  signalStatus(optimStatus(opt), errors.only=TRUE)
  
  ## .......... extract fitted parameters .....................
  if(!usecanonical) {
    optpar.canon <- NULL
    optpar.human <- opt$par
    names(optpar.human) <- names(startpar.human)
  } else {
    optpar.canon <- opt$par
    names(optpar.canon) <- names(startpar)
    optpar.human <- tohuman(optpar.canon)
    names(optpar.human) <- names(startpar.human)
  }
  opt$par       <- optpar.human
  opt$par.canon <- optpar.canon

  ## Finish in DPP case
  if(!is.null(DPP)){
    # all info that depends on the fitting method:
    Fit <- list(method    = "clik2",
                clfit     = opt,
                weightfun = weightfun,
                rmax      = rmax,
                objfun    = obj,
                objargs   = objargs,
                maxlogcl  = opt$value)
    # pack up
    clusters <- update(clusters, as.list(opt$par))
    result <- list(Xname      = Xname,
                   X          = X,
                   stationary = stationary,
                   fitted     = clusters,
                   modelname  = modelname,
                   po         = po,
                   lambda     = lambda,
                   Fit        = Fit)
    if(SAVE) attr(result, "h") <- get("h", envir=saveplace)
    return(result)
  }
  
  ## meaningful model parameters
  modelpar <- info$interpret(optpar.human, lambda)
  # infer parameter 'mu'
  if(isPCP) {
    # Poisson cluster process: extract parent intensity kappa
    kappa <- optpar.human[["kappa"]]
    # mu = mean cluster size
    mu <- if(stationary) lambda/kappa else eval.im(lambda/kappa)
  } else {
    # LGCP: extract variance parameter sigma2
    sigma2 <- optpar.human[["sigma2"]]
    # mu = mean of log intensity 
    mu <- if(stationary) log(lambda) - sigma2/2 else
          eval.im(log(lambda) - sigma2/2)    
  }
  # all info that depends on the fitting method:
  Fit <- list(method    = "clik2",
              clfit     = opt,
              weightfun = weightfun,
              rmax      = rmax,
              objfun    = obj,
              objargs   = objargs,
              maxlogcl  = opt$value)
  # pack up
  result <- list(Xname      = Xname,
                 X          = X,
                 stationary = stationary,
                 clusters   = clusters,
                 modelname  = modelname,
                 isPCP      = isPCP,
                 po         = po,
                 lambda     = lambda,
                 mu         = mu,
                 par        = optpar.human,
                 par.canon  = optpar.canon,
                 clustpar   = info$checkpar(par=optpar.human, old=FALSE, strict=strict),
                 clustargs  = info$checkclustargs(clargs$margs, old=FALSE, strict=strict), #clargs$margs,
                 modelpar   = modelpar,
                 covmodel   = clargs,
                 Fit        = Fit)
  if(SAVE) attr(result, "h") <- get("h", envir=saveplace)
  return(result)
}


kppmPalmLik <- function(X, Xname, po, clusters, control, weightfun, rmax,
                        algorithm="Nelder-Mead", DPP=NULL, ..., pint=NULL) {
  W <- as.owin(X)
  if(is.null(rmax))
    rmax <- rmax.rule("K", W, intensity(X))
  # identify pairs of points that contribute
  cl <- closepairs(X, rmax)
#  I <- cl$i
  J <- cl$j
  dIJ <- cl$d
  # compute weights for pairs of points
  if(is.function(weightfun)) {
    wIJ <- weightfun(dIJ)
#    sumweight <- sum(wIJ)
  } else {
    npairs <- length(dIJ)
    wIJ <- rep.int(1, npairs)
#    sumweight <- npairs
  }
  # convert window to mask, saving other arguments for later
  dcm <- do.call.matched(as.mask,
                         append(list(w=W), list(...)),
                         sieve=TRUE)
  M         <- dcm$result
  otherargs <- dcm$otherargs

  ## Detect DPP usage
  isDPP <- inherits(clusters, "detpointprocfamily")

  # compute intensity at data points
  # and c.d.f. of interpoint distance in window
  if(stationary <- is.stationary(po)) {
    # stationary unmarked Poisson process
    lambda <- intensity(X)
    lambdaJ <- rep(lambda, length(J))
    # compute cdf of distance between a uniform random point in W
    # and a randomly-selected point in X 
    g <- distcdf(X, M)
    # scaling constant is (integral of intensity) * (number of points)
    gscale <- npoints(X)^2
  } else {
    # compute fitted intensity at data points and in window
    lambdaX <- fitted(po, dataonly=TRUE)
    lambda <- lambdaM <- predict(po, locations=M)
    lambdaJ <- lambdaX[J] 
    # compute cdf of distance between a uniform random point in X 
    # and a random point in W with density proportional to intensity function
    g <- distcdf(X, M, dV=lambdaM)
    # scaling constant is (integral of intensity) * (number of points)
    gscale <- safevalue(integral.im(lambdaM) * npoints(X),
                        default=npoints(X)^2)
  }

  # Detect DPP model and change clusters and intensity correspondingly
  isDPP <- !is.null(DPP)
  if(isDPP){
    tmp <- dppmFixIntensity(DPP, lambda, po)
    clusters <- tmp$clusters
    lambda <- tmp$lambda
    po <- tmp$po
  }

  # trim 'g' to [0, rmax] 
  g <- g[with(g, .x) <= rmax,]
  # get pair correlation function (etc) for model
  info <- spatstatClusterModelInfo(clusters)
  pcfun      <- info$pcf
  funaux     <- info$funaux
  selfstart  <- info$selfstart
  isPCP      <- info$isPCP
  parhandler <- info$parhandler
  modelname  <- info$modelname
  # Assemble information required for computing pair correlation
  pcfunargs <- list(funaux=funaux)
  if(is.function(parhandler)) {
    # Additional parameters of cluster model are required.
    # These may be given as individual arguments,
    # or in a list called 'covmodel'
    clustargs <- if("covmodel" %in% names(otherargs))
                 otherargs[["covmodel"]] else otherargs
    clargs <- do.call(parhandler, clustargs)
    pcfunargs <- append(clargs, pcfunargs)
  } else clargs <- NULL
  # determine starting parameter values
  startpar <- selfstart(X)
  #' ............ experimental .........................
  strict <- !isFALSE(spatstat.options("kppm.strict"))
  if(!strict) pcfunargs <- append(pcfunargs, list(strict=strict))
  #' ............ experimental .........................
  usecanonical <- spatstat.options("kppm.canonical")
  if(usecanonical) {
     tocanonical <- info$tocanonical
     tohuman <- info$tohuman
     if(is.null(tocanonical) || is.null(tohuman)) {
       warning("Canonical parameters are not yet supported for this model")
       usecanonical <- FALSE
     }
  }
  startpar.human <- startpar
  if(usecanonical) {
    pcftheo <- pcfun
    startpar <- tocanonical(startpar)
    pcfun <- function(par, ...) { pcftheo(tohuman(par), ...) }
  }
  #' ............ experimental/debugger .........................
  whiu <- pint$whiu
  if(is.function(whiu) && usecanonical) {
    whiu.human <- whiu
    whiu <- function(par, ...) { whiu.human(tohuman(par), ...) }
  }
  TRACE <- isTRUE(pint$trace)
  if(SAVE <- isTRUE(pint$save)) {
    saveplace <- new.env()
    assign("h", NULL, envir=saveplace)
  } else saveplace <- NULL
  # .....................................................
  # create local function to evaluate pair correlation
  #  (with additional parameters 'pcfunargs' in its environment)
  paco <- function(d, par) {
    do.call(pcfun, append(list(par=par, rvals=d), pcfunargs))
  }
  # define objective function 
  if(!is.function(weightfun)) {
    # pack up necessary information
    objargs <- list(dIJ=dIJ, g=g, gscale=gscale,
                    sumloglam=safevalue(sum(log(lambdaJ))),
                    envir=environment(paco),
                    whiu=NULL,   # updated below
                    TRACE=FALSE, # updated below
                    saveplace=NULL, # updated below
                    BIGVALUE=1, # updated below
                    SMALLVALUE=.Machine$double.eps)
    # define objective function (with 'paco' in its environment)
    # This is the log Palm likelihood
    obj <- function(par, objargs) {
      with(objargs, {
        integ <- unlist(stieltjes(paco, g, par=par))
        integ <- pmax(SMALLVALUE, integ)
        logplik <- safevalue(sumloglam + sum(log(safevalue(paco(dIJ, par))))
                           - gscale * integ,
                           default=-BIGVALUE)
        if(is.function(whiu)) logplik <- whiu(par, logplik, 1)
        ## debugger
        if(isTRUE(TRACE)) {
          cat("Parameters:", fill=TRUE)
          print(par)
          splat("log Palm likelihood:", logplik)
        }
        if(is.environment(saveplace)) {
          h <- get("h", envir=saveplace)
          hplus <- as.data.frame(append(par, list(logplik=logplik)))
          h <- rbind(h, hplus)
          assign("h", h, envir=saveplace)
        }
        return(logplik)
      },
      enclos=objargs$envir)
    }
    ## determine a suitable large number to replace Inf
    objargs$BIGVALUE <- bigvaluerule(obj, objargs, startpar)
    objargs$whiu <- whiu
    objargs$saveplace <- saveplace
    objargs$TRACE <- TRACE
  } else {
    # create local function to evaluate  pair correlation(d) * weight(d)
    #  (with additional parameters 'pcfunargs', 'weightfun' in its environment)
    force(weightfun)
    wpaco <- function(d, par) {
      y <- do.call(pcfun, append(list(par=par, rvals=d), pcfunargs))
      w <- weightfun(d)
      return(y * w)
    }
    # pack up necessary information
    objargs <- list(dIJ=dIJ, wIJ=wIJ, g=g, gscale=gscale,
                    wsumloglam=safevalue(sum(wIJ * safevalue(log(lambdaJ)))),
                    envir=environment(wpaco),
                    whiu=NULL,   # updated below
                    TRACE=FALSE, # updated below
                    saveplace=NULL, # updated below
                    BIGVALUE=1, # updated below
                    SMALLVALUE=.Machine$double.eps)
    # define objective function (with 'paco', 'wpaco' in its environment)
    # This is the log Palm likelihood
    obj <- function(par, objargs) {
      with(objargs, {
        integ <- unlist(stieltjes(wpaco, g, par=par))
        integ <- pmax(SMALLVALUE, integ)
        logplik <- safevalue(wsumloglam +
                             sum(wIJ * log(safevalue(paco(dIJ, par))))
                             - gscale * integ,
                             default=-BIGVALUE)
        ## debugger
        if(isTRUE(TRACE)) {
          cat("Parameters:", fill=TRUE)
          print(par)
          splat("log Palm likelihood:", logplik)
        }
        if(is.environment(saveplace)) {
          h <- get("h", envir=saveplace)
          hplus <- as.data.frame(append(par, list(logplik=logplik)))
          h <- rbind(h, hplus)
          assign("h", h, envir=saveplace)
        }
        return(logplik)
        
      },
      enclos=objargs$envir)
    }
    ## determine a suitable large number to replace Inf
    objargs$BIGVALUE <- bigvaluerule(obj, objargs, startpar)
    objargs$whiu <- whiu
    objargs$saveplace <- saveplace
    objargs$TRACE <- TRACE
  }
  # arguments for optimization
  ctrl <- resolve.defaults(list(fnscale=-1), control, list(trace=0))
  optargs <- list(par=startpar, fn=obj, objargs=objargs,
                  control=ctrl, method=algorithm)
  ## DPP resolving algorithm and checking startpar
  changealgorithm <- length(startpar)==1 && algorithm=="Nelder-Mead"
  if(isDPP){
    alg <- dppmFixAlgorithm(algorithm, changealgorithm, clusters,
                            startpar.human)
    algorithm <- optargs$method <- alg$algorithm
    if(algorithm=="Brent" && changealgorithm){
      optargs$lower <- alg$lower
      optargs$upper <- alg$upper
    }
  }
  # optimize it
  opt <- do.call(optim, optargs)
  # raise warning/error if something went wrong
  signalStatus(optimStatus(opt), errors.only=TRUE)
  # Extract optimal values of parameters
  if(!usecanonical) {
    optpar.canon <- NULL
    optpar.human <- opt$par
    names(optpar.human) <- names(startpar.human)
  } else {
    optpar.canon <- opt$par
    names(optpar.canon) <- names(startpar)
    optpar.human <- tohuman(optpar.canon)
    names(optpar.human) <- names(startpar.human)
  }
  # Finish in DPP case
  if(!is.null(DPP)){
    opt$par <- optpar.human
    opt$par.canon <- optpar.canon
    # all info that depends on the fitting method:
    Fit <- list(method    = "palm",
                clfit     = opt,
                weightfun = weightfun,
                rmax      = rmax,
                objfun    = obj,
                objargs   = objargs,
                maxlogcl  = opt$value)
    # pack up
    clusters <- update(clusters, as.list(optpar.human))
    result <- list(Xname      = Xname,
                   X          = X,
                   stationary = stationary,
                   fitted     = clusters,
                   modelname  = modelname,
                   po         = po,
                   lambda     = lambda,
                   Fit        = Fit)
    if(SAVE) attr(result, "h") <- get("h", envir=saveplace)
    return(result)
  }
  # meaningful model parameters
  modelpar <- info$interpret(optpar.human, lambda)
  # infer parameter 'mu'
  if(isPCP) {
    # Poisson cluster process: extract parent intensity kappa
    kappa <- optpar.human[["kappa"]]
    # mu = mean cluster size
    mu <- if(stationary) lambda/kappa else eval.im(lambda/kappa)
  } else {
    # LGCP: extract variance parameter sigma2
    sigma2 <- optpar.human[["sigma2"]]
    # mu = mean of log intensity 
    mu <- if(stationary) log(lambda) - sigma2/2 else
          eval.im(log(lambda) - sigma2/2)    
  }
  # all info that depends on the fitting method:
  Fit <- list(method    = "palm",
              clfit     = opt,
              weightfun = weightfun,
              rmax      = rmax,
              objfun    = obj,
              objargs   = objargs,
              maxlogcl  = opt$value)
  # pack up
  result <- list(Xname      = Xname,
                 X          = X,
                 stationary = stationary,
                 clusters   = clusters,
                 modelname  = modelname,
                 isPCP      = isPCP,
                 po         = po,
                 lambda     = lambda,
                 mu         = mu,
                 par        = optpar.human,
                 par.canon  = optpar.canon,
                 clustpar   = info$checkpar(par=optpar.human, old=FALSE, strict=strict),
                 clustargs  = info$checkclustargs(clargs$margs, old=FALSE, strict=strict), #clargs$margs,
                 modelpar   = modelpar,
                 covmodel   = clargs,
                 Fit        = Fit)
  if(SAVE) attr(result, "h") <- get("h", envir=saveplace)
  return(result)
}

improve.kppm <- local({

  fnc <- function(r, eps, g){ (g(r) - 1)/(g(0) - 1) - eps}

  improve.kppm <- function(object, type=c("quasi", "wclik1", "clik1"),
                           rmax = NULL, eps.rmax = 0.01,
                           dimyx = 50, maxIter = 100, tolerance = 1e-06,
                           fast = TRUE, vcov = FALSE, fast.vcov = FALSE,
                           verbose = FALSE,
                           save.internals = FALSE) {
    verifyclass(object, "kppm")
    type <- match.arg(type)
    gfun <- pcfmodel(object)
    X <- object$X
    win <- as.owin(X)
    ## simple (rectangular) grid quadrature scheme
    ## (using pixels with centers inside owin only)
    mask <- as.mask(win, dimyx = dimyx)
    wt <- pixellate(win, W = mask)
    wt <- wt[mask]
    Uxy <- rasterxy.mask(mask)
    U <- ppp(Uxy$x, Uxy$y, window = win, check=FALSE)
    U <- U[mask]
#    nU <- npoints(U)
    Yu <- pixellate(X, W = mask)
    Yu <- Yu[mask]
    
    ## covariates at quadrature points
    po <- object$po
    Z <- model.images(po, mask)
    Z <- sapply(Z, "[", i=U)

    ##obtain initial beta estimate using composite likelihood
    beta0 <- coef(po)
    
    ## determining the dependence range
    if (type != "clik1" && is.null(rmax))
      {
        diamwin <- diameter(win)
        rmax <- if(fnc(diamwin, eps.rmax, gfun) >= 0) diamwin else
                uniroot(fnc, lower = 0, upper = diameter(win),
                        eps=eps.rmax, g=gfun)$root
        if(verbose)
          splat(paste0("type: ", type, ", ",
                       "dependence range: ", rmax, ", ",
                       "dimyx: ", dimyx, ", g(0) - 1:", gfun(0) -1))
      }
    ## preparing the WCL case
    if (type == "wclik1")
      Kmax <- 2*pi * integrate(function(r){r * (gfun(r) - 1)},
                               lower=0, upper=rmax)$value * exp(c(Z %*% beta0))
    ## the g()-1 matrix without tapering
    if (!fast || (vcov && !fast.vcov)){
      if (verbose)
        cat("computing the g(u_i,u_j)-1 matrix ...")
      gminus1 <- matrix(gfun(c(pairdist(U))) - 1, U$n, U$n)
      if (verbose)
        cat("..Done.\n")
    }
    
    if ( (fast && type == "quasi") | fast.vcov ){
      if (verbose)
        cat("computing the sparse G-1 matrix ...\n")
      ## Non-zero gminus1 entries (when using tapering)
      cp <- crosspairs(U,U,rmax,what="ijd")
      if (verbose)
        cat("crosspairs done\n")
      Gtap <- (gfun(cp$d) - 1)
      if(vcov){
        if(fast.vcov){
          gminus1 <- Matrix::sparseMatrix(i=cp$i, j=cp$j,
                                          x=Gtap, dims=c(U$n, U$n))
        } else{
          if(fast)
            gminus1 <- matrix(gfun(c(pairdist(U))) - 1, U$n, U$n)
        }
      }
      if (verbose & type!="quasi")
        cat("..Done.\n")
    }
       
    if (type == "quasi" && fast){
      mu0 <- exp(c(Z %*% beta0)) * wt
      mu0root <- sqrt(mu0)
      sparseG <- Matrix::sparseMatrix(i=cp$i, j=cp$j,
                                      x=mu0root[cp$i] * mu0root[cp$j] * Gtap,
                                      dims=c(U$n, U$n))
      Rroot <- Matrix::Cholesky(sparseG, perm = TRUE, Imult = 1)
      ##Imult=1 means that we add 1*I
      if (verbose)
        cat("..Done.\n")
    }
    
    ## iterative weighted least squares/Fisher scoring
    bt <- beta0
    noItr <- 1
    repeat {
      mu <- exp(c(Z %*% bt)) * wt
      mu.root <- sqrt(mu)
      ## the core of estimating equation: ff=phi
      ## in case of ql, \phi=V^{-1}D=V_\mu^{-1/2}x where (G+I)x=V_\mu^{1/2} Z
      ff <- switch(type,
                   clik1 = Z,
                   wclik1= Z/(1 + Kmax),
                   quasi = if(fast){
                     Matrix::solve(Rroot, mu.root * Z)/mu.root
                   } else{
                     solve(diag(U$n) + t(gminus1 * mu), Z)
                   }
                   )
      ##alternative
      ##R=chol(sparseG+sparseMatrix(i=c(1:U$n),j=c(1:U$n),
      ##                            x=rep(1,U$n),dims=c(U$n,U$n)))
      ##ff2 <- switch(type,
      ##              clik1 = Z,
      ##              wclik1= Z/(1 + Kmax),
      ##              quasi = if (fast)
      ##                         solve(R,solve(t(R), mu.root * Z))/mu.root
      ##                      else solve(diag(U$n) + t(gminus1 * mu), Z))
      ## print(summary(as.numeric(ff)-as.numeric(ff2)))
      ## the estimating equation: u_f(\beta)
      uf <- (Yu - mu) %*% ff
      ## inverse of minus expectation of Jacobian matrix: I_f
      Jinv <- solve(t(Z * mu) %*% ff)
      if(maxIter==0){
        ## This is a built-in early exit for vcov internal calculations
        break
      }
      deltabt <- as.numeric(uf %*% Jinv)
      if (any(!is.finite(deltabt))) {
        warning(paste("Infinite value, NA or NaN appeared",
                      "in the iterative weighted least squares algorithm.",
                      "Returning the initial intensity estimate unchanged."),
                call.=FALSE)
        return(object)
      }
      ## updating the present estimate of \beta
      bt <- bt + deltabt
      if (verbose)
        splat(paste0("itr: ", noItr, ",\nu_f: ", as.numeric(uf),
                     "\nbeta:", bt, "\ndeltabeta:", deltabt))
      if (max(abs(deltabt/bt)) <= tolerance || max(abs(uf)) <= tolerance)
        break
      if (noItr > maxIter)
        stop("Maximum number of iterations reached without convergence.")
      noItr <- noItr + 1
    }
    out <- object
    out$po$coef.orig <- beta0
    out$po$coef <- bt
    loc <- if(is.sob(out$lambda)) as.mask(out$lambda) else mask
    out$lambda <- predict(out$po, locations = loc)
    out$improve <- list(type = type,
                        rmax = rmax,
                        dimyx = dimyx,
                        fast = fast,
                        fast.vcov = fast.vcov)
    if(save.internals){
      out$improve <- append(out$improve, list(ff=ff, uf=uf, J.inv=Jinv))
    }
    if(vcov){
      if (verbose)
        cat("computing the asymptotic variance ...\n")
      ## variance of the estimation equation: Sigma_f = Var(u_f(bt))
      trans <- if(fast) Matrix::t else t
      Sig <- trans(ff) %*% (ff * mu) + trans(ff * mu) %*% gminus1 %*% (ff * mu)
      ## note Abdollah's G does not have mu.root inside...
      ## the asymptotic variance of \beta:
      ##         inverse of the Godambe information matrix
      out$vcov <- as.matrix(Jinv %*% Sig %*% Jinv)
    }
    return(out)
  }
  improve.kppm
})


is.kppm <- function(x) { inherits(x, "kppm")}

print.kppm <- print.dppm <- function(x, ...) {

  isPCP <- x$isPCP
  # detect DPP
  isDPP <- inherits(x, "dppm")
  # handle outdated objects - which were all cluster processes
  if(!isDPP && is.null(isPCP)) isPCP <- TRUE

  terselevel <- spatstat.options('terse')
  digits <- getOption('digits')
  
  splat(if(x$stationary) "Stationary" else "Inhomogeneous",
        if(isDPP) "determinantal" else if(isPCP) "cluster" else "Cox",
        "point process model")

  Xname <- x$Xname
  if(waxlyrical('extras', terselevel) && nchar(Xname) < 20) {
    has.subset <- ("subset" %in% names(x$call))
    splat("Fitted to",
          if(has.subset) "(a subset of)" else NULL,
          "point pattern dataset", sQuote(Xname))
  }

  if(waxlyrical('gory', terselevel)) {
    switch(x$Fit$method,
           mincon = {
             splat("Fitted by minimum contrast")
             splat("\tSummary statistic:", x$Fit$StatName)
           },
           clik  =,
           clik2 = {
             splat("Fitted by maximum second order composite likelihood")
             splat("\trmax =", x$Fit$rmax)
             if(!is.null(wtf <- x$Fit$weightfun)) {
               a <- attr(wtf, "selfprint") %orifnull% pasteFormula(wtf)
               splat("\tweight function:", a)
             }
           },
           palm = {
             splat("Fitted by maximum Palm likelihood")
             splat("\trmax =", x$Fit$rmax)
             if(!is.null(wtf <- x$Fit$weightfun)) {
               a <- attr(wtf, "selfprint") %orifnull% pasteFormula(wtf)
               splat("\tweight function:", a)
             }
           },
           warning(paste("Unrecognised fitting method", sQuote(x$Fit$method)))
           )
  }

  parbreak(terselevel)
  
  # ............... trend .........................

  if(!(isDPP && is.null(x$fitted$intensity)))
    print(x$po, what="trend")

  # ..................... clusters ................

  # DPP case
  if(isDPP){
      splat("Fitted DPP model:")
      print(x$fitted)
      return(invisible(NULL))
  }

  tableentry <- spatstatClusterModelInfo(x$clusters)
  
  splat(if(isPCP) "Cluster" else "Cox",
        "model:", tableentry$printmodelname(x))
  cm <- x$covmodel
  if(!isPCP) {
    # Covariance model - LGCP only
    splat("\tCovariance model:", cm$model)
    margs <- cm$margs
    if(!is.null(margs)) {
      nama <- names(margs)
      tags <- ifelse(nzchar(nama), paste(nama, "="), "")
      tagvalue <- paste(tags, margs)
      splat("\tCovariance parameters:",
            paste(tagvalue, collapse=", "))
    }
  }
  pc <- x$par.canon
  if(!is.null(pc)) {
    splat("Fitted canonical parameters:")
    print(pc, digits=digits)
  }
  pa <- x$clustpar
  if (!is.null(pa)) {
    splat("Fitted",
          if(isPCP) "cluster" else "covariance",
          "parameters:")
    print(pa, digits=digits)
  }

  if(!is.null(mu <- x$mu)) {
    if(isPCP) {
      splat("Mean cluster size: ",
            if(!is.im(mu)) paste(signif(mu, digits), "points") else "[pixel image]")
    } else {
      splat("Fitted mean of log of random intensity:",
            if(!is.im(mu)) signif(mu, digits) else "[pixel image]")
    }
  }
  if(isDPP) {
    rx <- repul(x)
    splat(if(is.im(rx)) "(Average) strength" else "Strength",
          "of repulsion:", signif(mean(rx), 4))
  }
  invisible(NULL)
}

plot.kppm <- local({

  plotem <- function(x, ..., main=dmain, dmain) { plot(x, ..., main=main) }
  
  plot.kppm <- function(x, ...,
                        what=c("intensity", "statistic", "cluster"),
                        pause=interactive(),
                        xname) {
    ## catch objectname from dots if present otherwise deparse x:
    if(missing(xname)) xname <- short.deparse(substitute(x))
    nochoice <- missing(what)
    what <- pickoption("plot type", what,
                       c(statistic="statistic",
                         intensity="intensity",
                         cluster="cluster"),
                       multi=TRUE)
    ## handle older objects
    Fit <- x$Fit
    if(is.null(Fit)) {
      warning("kppm object is in outdated format")
      Fit <- x
      Fit$method <- "mincon"
    }
    ## Catch locations for clusters if given
    loc <- list(...)$locations
    inappropriate <- (nochoice & ((what == "intensity") & (x$stationary))) |
             ((what == "statistic") & (Fit$method != "mincon")) |
             ((what == "cluster") & (identical(x$isPCP, FALSE))) | 
             ((what == "cluster") & (!x$stationary) & is.null(loc))

    if(!nochoice && !x$stationary && "cluster" %in% what && is.null(loc))
      stop("Please specify additional argument ", sQuote("locations"),
           " which will be passed to the function ",
           sQuote("clusterfield"), ".")

    if(any(inappropriate)) {
      what <- what[!inappropriate]
      if(length(what) == 0){
        message("Nothing meaningful to plot. Exiting...")
        return(invisible(NULL))
      }
    }
    pause <- pause && (length(what) > 1)
    if(pause) opa <- par(ask=TRUE)
    for(style in what)
      switch(style,
             intensity={
               plotem(x$po, ...,
                      dmain=c(xname, "Intensity"),
                      how="image", se=FALSE)
             },
             statistic={
               plotem(Fit$mcfit, ...,
                      dmain=c(xname, Fit$StatName))
             },
             cluster={
               plotem(clusterfield(x, locations = loc, verbose=FALSE), ...,
                      dmain=c(xname, "Fitted cluster"))
             })
    if(pause) par(opa)
    return(invisible(NULL))
  }

  plot.kppm
})


predict.kppm <- predict.dppm <- function(object, ...) {
  se <- resolve.1.default(list(se=FALSE), list(...))
  interval <- resolve.1.default(list(interval="none"), list(...))
  if(se) warning("Standard error calculation assumes a Poisson process")
  if(interval != "none")
    warning(paste(interval, "interval calculation assumes a Poisson process"))
  predict(as.ppm(object), ...)
}

fitted.kppm <- fitted.dppm <- function(object, ...) {
  fitted(as.ppm(object), ...)
}

residuals.kppm <- residuals.dppm <- function(object, ...) {
  type <- resolve.1.default(list(type="raw"), list(...))
  if(type != "raw")
    warning(paste("calculation of", type, "residuals",
                  "assumes a Poisson process"))
  residuals(as.ppm(object), ...)
}

simulate.kppm <- function(object, nsim=1, seed=NULL, ...,
                          window=NULL, covariates=NULL,
                          verbose=TRUE, retry=10,
                          drop=FALSE) {
  starttime <- proc.time()
  verbose <- verbose && (nsim > 1)
  check.1.real(retry)
  # .... copied from simulate.lm ....
  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
    runif(1)
  if (is.null(seed))
    RNGstate <- get(".Random.seed", envir = .GlobalEnv)
  else {
    R.seed <- get(".Random.seed", envir = .GlobalEnv)
    set.seed(seed)
    RNGstate <- structure(seed, kind = as.list(RNGkind()))
    on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
  }
  
  # ..................................
  # determine window for simulation results
  if(!is.null(window)) {
    stopifnot(is.owin(window))
    win <- window
  } else {
    win <- as.owin(object)
  }
  # ..................................
  # determine parameters
  mp <- as.list(object$modelpar)

  # parameter 'mu'
  # = parent intensity of cluster process
  # = mean log intensity of log-Gaussian Cox process
  
  if(is.null(covariates) && (object$stationary || is.null(window))) {
    # use existing 'mu' (scalar or image)
    mu <- object$mu
  } else {
    # recompute 'mu' using new data
    switch(object$clusters,
           Cauchy=,
           VarGamma=,
           Thomas=,
           MatClust={
             # Poisson cluster process
             kappa <- mp$kappa
             lambda <- predict(object, window=win, covariates=covariates)
             mu <- eval.im(lambda/kappa)
           },
           LGCP={
             # log-Gaussian Cox process
             sigma2 <- mp$sigma2
             lambda <- predict(object, window=win, covariates=covariates)
             mu <- eval.im(log(lambda) - sigma2/2)
           },
           stop(paste("Simulation of", sQuote(object$clusters),
                      "processes is not yet implemented"))
           )
  }

  # prepare data for execution
  out <- list()
  switch(object$clusters,
         Thomas={
           kappa <- mp$kappa
           sigma <- mp$sigma
           cmd <- expression(rThomas(kappa,sigma,mu,win, ...))
           dont.complain.about(kappa, sigma, mu)
         },
         MatClust={
           kappa <- mp$kappa
           r     <- mp$R
           cmd   <- expression(rMatClust(kappa,r,mu,win, ...))
           dont.complain.about(kappa, r)
         },
         Cauchy = {
           kappa <- mp$kappa
           omega <- mp$omega
           cmd   <- expression(rCauchy(kappa, omega, mu, win, ...))
           dont.complain.about(kappa, omega, mu)
         },
         VarGamma = {
           kappa  <- mp$kappa
           omega  <- mp$omega
           nu.ker <- object$covmodel$margs$nu.ker
           cmd    <- expression(rVarGamma(kappa, nu.ker, omega, mu, win, ...))
           dont.complain.about(kappa, nu.ker, omega, mu)
         },
         LGCP={
           sigma2 <- mp$sigma2
           alpha  <- mp$alpha
           cm <- object$covmodel
           model <- cm$model
           margs <- cm$margs
           param <- append(list(var=sigma2, scale=alpha), margs)
           #' 
           if(!is.im(mu)) {
             # model will be simulated in 'win'
             cmd <- expression(rLGCP(model=model, mu=mu, param=param,
                               ..., win=win))
             #' check that RandomFields package recognises parameter format
             rfmod <- try(rLGCP(model, mu=mu, param=param, win=win,
                              ..., modelonly=TRUE))
           } else {
             # model will be simulated in as.owin(mu), then change window
             cmd <- expression(rLGCP(model=model, mu=mu, param=param,
                               ...)[win])
             #' check that RandomFields package recognises parameter format
             rfmod <- try(rLGCP(model, mu=mu, param=param, 
                              ..., modelonly=TRUE))
           }
           #' suppress warnings from code checker
           dont.complain.about(model, mu, param)
           #' check that model is recognised
           if(inherits(rfmod, "try-error"))
             stop(paste("Internal error in simulate.kppm:",
                        "unable to build Random Fields model",
                        "for log-Gaussian Cox process"))
         })
  
  # run
  if(verbose) {
    cat(paste("Generating", nsim, "simulations... "))
    state <- list()
  }
  for(i in 1:nsim) {
    out[[i]] <- try(eval(cmd))
    if(verbose) state <- progressreport(i, nsim, state=state)
  }
  # detect failures
  if(any(bad <- unlist(lapply(out, inherits, what="try-error")))) {
    nbad <- sum(bad)
    gripe <- paste(nbad,
                   ngettext(nbad, "simulation was", "simulations were"),
                   "unsuccessful")
    if(verbose) splat(gripe)
    if(retry <= 0) {
      fate <- "returned as NULL"
      out[bad] <- list(NULL)
    } else {
      if(verbose) cat("Retrying...")
      ntried <- 0
      while(ntried < retry) {
        ntried <- ntried + 1
        for(j in which(bad))
          out[[j]] <- try(eval(cmd))
        bad <- unlist(lapply(out, inherits, what="try-error"))
        nbad <- sum(bad)
        if(nbad == 0) break
      }
      if(verbose) cat("Done.\n")
      fate <- if(nbad == 0) "all recomputed" else
              paste(nbad, "simulations still unsuccessful")
      fate <- paste(fate, "after", ntried,
                    ngettext(ntried, "further try", "further tries"))
    }
    warning(paste(gripe, fate, sep=": "))
  }
  if(verbose)
    cat("Done.\n")
  #' pack up
  out <- simulationresult(out, nsim, drop)
  out <- timed(out, starttime=starttime)
  attr(out, "seed") <- RNGstate
  return(out)
}

formula.kppm <- formula.dppm <- function(x, ...) {
  formula(x$po, ...)
}

terms.kppm <- terms.dppm <- function(x, ...) {
  terms(x$po, ...)
}

labels.kppm <- labels.dppm <- function(object, ...) {
  labels(object$po, ...)
}

update.kppm <- function(object, ..., evaluate=TRUE) {
  argh <- list(...)
  nama <- names(argh)
  callframe <- object$callframe
  envir <- environment(terms(object))
  #' look for a formula argument
  fmla <- formula(object)
  jf <- integer(0)
  if(!is.null(trend <- argh$trend)) {
    if(!can.be.formula(trend))
      stop("Argument \"trend\" should be a formula")
    fmla <- newformula(formula(object), trend, callframe, envir)
    jf <- which(nama == "trend")
  } else if(any(isfo <- sapply(argh, can.be.formula))) {
    if(sum(isfo) > 1) {
      if(!is.null(nama)) isfo <- isfo & nzchar(nama)
      if(sum(isfo) > 1)
        stop(paste("Arguments not understood:",
                   "there are two unnamed formula arguments"))
    }
    jf <- which(isfo)
    fmla <- argh[[jf]]
    fmla <- newformula(formula(object), fmla, callframe, envir)
  }

  #' look for a point pattern or quadscheme
  if(!is.null(X <- argh$X)) {
    if(!inherits(X, c("ppp", "quad")))
      stop(paste("Argument X should be a formula,",
                 "a point pattern or a quadrature scheme"))
    jX <- which(nama == "X")
  } else if(any(ispp <- sapply(argh, inherits, what=c("ppp", "quad")))) {
    if(sum(ispp) > 1) {
      if(!is.null(nama)) ispp <- ispp & nzchar(nama)
      if(sum(ispp) > 1)
        stop(paste("Arguments not understood:",
                   "there are two unnamed point pattern/quadscheme arguments"))
    }
    jX <- which(ispp)
    X <- argh[[jX]]
  } else {
    X <- object$X
    jX <- integer(0)
  }
  Xexpr <- if(length(jX) > 0) sys.call()[[2L + jX]] else NULL
  #' remove arguments just recognised, if any
  jused <- c(jf, jX)
  if(length(jused) > 0) {
    argh <- argh[-jused]
    nama <- names(argh)
  }
  #' update the matched call
  thecall <- getCall(object)
  methodname <- as.character(thecall[[1L]])
  switch(methodname,
         kppm.formula = {
	   # original call has X = [formula with lhs]
	   if(!is.null(Xexpr)) {
	     lhs.of.formula(fmla) <- Xexpr
	   } else if(is.null(lhs.of.formula(fmla))) {
	     lhs.of.formula(fmla) <- as.name('.')
	   }
           oldformula <- as.formula(getCall(object)$X)
           thecall$X <- newformula(oldformula, fmla, callframe, envir)
         },
         {
	   # original call has X = ppp and trend = [formula without lhs]
           oldformula <- as.formula(getCall(object)$trend %orifnull% (~1))
	   fom <-  newformula(oldformula, fmla, callframe, envir)
	   if(!is.null(Xexpr))
	      lhs.of.formula(fom) <- Xexpr
	   if(is.null(lhs.of.formula(fom))) {
	      # new call has same format
	      thecall$trend <- fom
  	      if(length(jX) > 0)
  	        thecall$X <- X
	   } else {
	      # new call has formula with lhs
	      thecall$trend <- NULL
	      thecall$X <- fom
	   }
         })
  knownnames <- unique(c(names(formals(kppm.ppp)),
                         names(formals(mincontrast)),
                         names(formals(optim))))
  knownnames <- setdiff(knownnames,
                        c("X", "trend",
                          "observed", "theoretical",
                          "fn", "gr", "..."))
  ok <- nama %in% knownnames
  thecall <- replace(thecall, nama[ok], argh[ok])
  thecall$formula <- NULL # artefact of 'step', etc
  thecall[[1L]] <- as.name("kppm")
  if(!evaluate)
    return(thecall)
  out <- eval(thecall, envir=parent.frame(), enclos=envir)
  #' update name of data
  if(length(jX) == 1) {
    mc <- match.call()
    Xlang <- mc[[2L+jX]]
    out$Xname <- short.deparse(Xlang)
  }
  #'
  return(out)
}

unitname.kppm <- unitname.dppm <- function(x) {
  return(unitname(x$X))
}

"unitname<-.kppm" <- "unitname<-.dppm" <- function(x, value) {
  unitname(x$X) <- value
  if(!is.null(x$Fit$mcfit)) {
    unitname(x$Fit$mcfit) <- value
  } else if(is.null(x$Fit)) {
    warning("kppm object in outdated format")
    if(!is.null(x$mcfit))
      unitname(x$mcfit) <- value
  }
  return(x)
}

as.fv.kppm <- as.fv.dppm <- function(x) {
  if(x$Fit$method == "mincon")
    return(as.fv(x$Fit$mcfit))
  gobs <- pcfinhom(x$X, lambda=x, correction="good", update=FALSE)
  gfit <- (pcfmodel(x))(gobs$r)
  g <- bind.fv(gobs,
               data.frame(fit=gfit), 
               "%s[fit](r)",
               "predicted %s for fitted model")
  return(g)
}

coef.kppm <- coef.dppm <- function(object, ...) {
  return(coef(object$po))
}


Kmodel.kppm <- function(model, ...) {
  Kpcf.kppm(model, what="K")
}

pcfmodel.kppm <- function(model, ...) {
  Kpcf.kppm(model, what="pcf")
}

Kpcf.kppm <- function(model, what=c("K", "pcf", "kernel")) {
  what <- match.arg(what)
  # Extract function definition from internal table
  clusters <- model$clusters
  tableentry <- spatstatClusterModelInfo(clusters)
  if(is.null(tableentry))
    stop("No information available for", sQuote(clusters), "cluster model")
  fun <- tableentry[[what]]
  if(is.null(fun))
    stop("No expression available for", what, "for", sQuote(clusters),
         "cluster model")
  # Extract model parameters
  par <- model$par
  # Extract auxiliary definitions (if applicable)
  funaux <- tableentry$funaux
  # Extract covariance model (if applicable)
  cm <- model$covmodel
  model <- cm$model
  margs <- cm$margs
  #
  f <- function(r) as.numeric(fun(par=par, rvals=r,
                                  funaux=funaux, model=model, margs=margs))
  return(f)
}

is.stationary.kppm <- is.stationary.dppm <- function(x) {
  return(x$stationary)
}

is.poisson.kppm <- function(x) {
  switch(x$clusters,
         Cauchy=,
         VarGamma=,
         Thomas=,
         MatClust={
           # Poisson cluster process
           mu <- x$mu
           return(!is.null(mu) && (max(mu) == 0))
         },
         LGCP = {
           # log-Gaussian Cox process
           sigma2 <- x$par[["sigma2"]]
           return(sigma2 == 0)
         },
         return(FALSE))
}

# extract ppm component

as.ppm.kppm <- as.ppm.dppm <- function(object) {
  object$po
}

# other methods that pass through to 'ppm'

as.owin.kppm <- as.owin.dppm <- function(W, ..., from=c("points", "covariates"), fatal=TRUE) {
  from <- match.arg(from)
  as.owin(as.ppm(W), ..., from=from, fatal=fatal)
}

domain.kppm <- Window.kppm <- domain.dppm <-
  Window.dppm <- function(X, ..., from=c("points", "covariates")) {
  from <- match.arg(from)
  as.owin(X, from=from)
}

model.images.kppm <-
  model.images.dppm <- function(object, W=as.owin(object), ...) {
  model.images(as.ppm(object), W=W, ...)
}

model.matrix.kppm <-
  model.matrix.dppm <- function(object,
                                data=model.frame(object, na.action=NULL), ...,
                                Q=NULL, 
                                keepNA=TRUE) {
  if(missing(data)) data <- NULL
  model.matrix(as.ppm(object), data=data, ..., Q=Q, keepNA=keepNA)
}

model.frame.kppm <- model.frame.dppm <- function(formula, ...) {
  model.frame(as.ppm(formula), ...)
}

logLik.kppm <- logLik.dppm <- function(object, ...) {
  cl <- object$Fit$maxlogcl
  if(is.null(cl))
    stop(paste("logLik is only available for kppm objects fitted with",
               "method='palm' or method='clik2'"),
         call.=FALSE)
  ll <- logLik(as.ppm(object)) # to inherit class and d.f.
  ll[] <- cl
  return(ll)
}
 
AIC.kppm <- AIC.dppm <- function(object, ..., k=2) {
  cl <- logLik(object)
  df <- attr(cl, "df")
  return(- 2 * as.numeric(cl) + k * df)
}

extractAIC.kppm <- extractAIC.dppm <- function (fit, scale = 0, k = 2, ...) {
  cl <- logLik(fit)
  edf <- attr(cl, "df")
  aic <- - 2 * as.numeric(cl) + k * edf
  return(c(edf, aic))
}

nobs.kppm <- nobs.dppm <- function(object, ...) { nobs(as.ppm(object)) }

psib <- function(object) UseMethod("psib")

psib.kppm <- function(object) {
  clus <- object$clusters
  info <- spatstatClusterModelInfo(clus)
  if(!info$isPCP) {
    warning("The model is not a cluster process")
    return(NA)
  }
  g <- pcfmodel(object)
  p <- 1 - 1/g(0)
  return(p)
}
